program test_solver
    use FChaseSolver
    implicit none;
    integer :: n;
    real, pointer :: a(:);
    real, pointer :: b(:);
    real, pointer :: c(:);
    real, pointer :: d(:);
    real, pointer :: x(:);
    real, pointer :: AMat(:, :)

    n = 4;
    allocate(a(n - 1));
    allocate(b(n));
    allocate(c(n - 1));
    allocate(d(n));
    allocate(x(n));
    a = [1, 1, 1];
    b = [1, 2, 3, 4]
    c = [1, 1, 1]
    d = [3, 8, 15, 19]

    call  fcs_idep(x, n, a, b, c, d);

    print *, x;

    allocate(AMat(n, n))

    AMat = reshape([1, 1, 0, 0, &
            1, 2, 1, 0, &
            0, 1, 3, 1, &
            0, 0, 1, 4], shape = [n, n]);

    print *, "The solution should be [1., 2., 3., 4.]";
    print *, matmul(AMat, x)

    deallocate(a);
    deallocate(b);
    deallocate(c);
    deallocate(d);
    deallocate(x);

end program test_solver